function res = gdatest(stego, cover, ntrain, ntest, loop, verbose)
if nargin < 6
    verbose = 1;
end
[s1, h1] = size(stego);
[s2, h2] = size(cover);
ntrain = ntrain / 2;
ntest = ntest / 2;
if (ntrain+ntest > s1) || (ntrain + ntest > s2)
    error('number of training and testing samples must less than .');
end
if (h1 ~= h2) || (s1 ~= s2)
    error('dimension of stego / cover does not match.');
end

testlbl = [ones(ntest,1);zeros(ntest,1)];
res = struct('MEAN_ACC', 0, 'acc', zeros(loop, 1), 'train_time', zeros(loop,1));
for i = 1:loop
    if verbose
        fprintf('round %d:\n', i);
    end
    p = randperm(s1);
    trainset = [stego(p(1:ntrain),:);cover(p(1:ntrain),:)];
    testset = [stego(p(ntrain+1:ntrain+ntest),:);cover(p(ntrain+1:ntrain+ntest),:)];
    
    %scale
    model = struct('sigma',5.5);
    model.meantra=mean(trainset);
    model.stdtra=std(trainset);
    trainset = scaleData(trainset, model.meantra, model.stdtra);
    testset = scaleData(testset, model.meantra, model.stdtra);

    tic;
    model.dataGDA = buildGDA(trainset, [ntrain, ntrain], model.sigma);
    res.train_time(i) = toc;
    
    %map to kernel space
    fprintf('.');
    vtrain = gdamapping(trainset, trainset, model.dataGDA, 1, model.sigma);
    fprintf('.');
    vtest = gdamapping(testset, trainset, model.dataGDA, 1, model.sigma);
    fprintf('.');
    fprintf('\n');

    model.vStego = mean(vtrain(1:ntrain))
    model.vCover = mean(vtrain(ntrain+1:2*ntrain))
    model.mean = (model.vStego + model.vCover)/2
    model.ab = (model.vStego >= model.vCover);
    
    val = (vtest - model.mean) >= 0;
    if ~model.ab
        val = ~val;
    end
    res.acc(i) = sum(val == testlbl) / ntest / 2;
    
    %res.tp = sum(abs(v1(1:ntest))<abs(v2(1:ntest)))/ntest;
    %res.tn = sum(abs(v1(ntest+1:end))>=abs(v2(ntest+1:end)))/ntest;

    if verbose
        fprintf('acc = %g\n',res.acc(i));
        fprintf('time = %g\n',res.train_time(i));
    end
end
res.MEAN_ACC = mean(res.acc);
if verbose
    disp(res)
end

end

function data = scaleData(data, m, s)
%m - mean
%s - std
p = size(data, 1);
data = (data - repmat(m, p, 1))./repmat(s, p, 1);
end